/*
 * volumegenerator.cpp
 *
 * Copyright 2012 The University of Auckland
 *
 *  Created on: 2/07/2012
 *      Author: cvan680
 *
 *
 * TODO:
 *
 * 	- More options: e.g. Scalar data, Vector data, Tensor data, multifield volume...
 *
 */
//#include "custommodules/ccvethesis/datastructures/ctensor.h"

#include "volumegenerator.h"
#include "voreen/core/datastructures/volume/volumefactory.h"
#include "voreen/core/datastructures/volume/volumeatomic.h"
#include "voreen/core/datastructures/volume/volumeram.h"
#include "voreen/core/datastructures/volume/volume.h"
#include "tgt/textureunit.h"
#include "tgt/filesystem.h"


#include <Eigen/Core>
#include <Eigen/Eigen>


using tgt::vec3;
using tgt::svec3;

namespace voreen {

VolumeGenerator::VolumeGenerator()
	: VolumeProcessor()
	, outport(Port::OUTPORT, "volumecollection", 0)
	, dimension_x_prop("dimension_x", "Dimension (x)", 128, 1, 512)
	, dimension_y_prop("dimension_y", "Dimension (y)", 128, 1, 512)
	, dimension_z_prop("dimension_z", "Dimension (z)", 128, 1, 512)
	, dimensions_cubed_prop("dimensions_cubed", "Cube volume?", true)
	, num_scalar_volumes_prop("num_scalar_volumes", "Scalar volumes", 0, 0, 6)
	, num_vector_volumes_prop("num_vector_volumes", "Vector volumes", 0, 0, 6)
	, num_tensor_volumes_prop("num_tensor_volumes", "Tensor volumes", 0, 0, 6)
	, update_prop("update", "Update outputs")
{
    addPort(outport);

    /* TODO:
     * - maybe offer several different generating functions..?
     * - should we support tensor datasets at all?
     * - where do we build up the volumes?
     */

	dimension_x_prop.onChange(CallMemberAction<VolumeGenerator>(this, &VolumeGenerator::onXDimensionsChange));
	dimension_y_prop.onChange(CallMemberAction<VolumeGenerator>(this, &VolumeGenerator::onYDimensionsChange));
	dimension_z_prop.onChange(CallMemberAction<VolumeGenerator>(this, &VolumeGenerator::onZDimensionsChange));

    addProperty(dimension_x_prop);
    addProperty(dimension_y_prop);
    addProperty(dimension_z_prop);
    addProperty(dimensions_cubed_prop);


    addProperty(num_scalar_volumes_prop);
    addProperty(num_vector_volumes_prop);
    addProperty(num_tensor_volumes_prop);

	update_prop.onClick(CallMemberAction<VolumeGenerator>(this, &VolumeGenerator::update));
	addProperty(update_prop);



/* TEMP! just checking out if it works
	using namespace Eigen;


	MatrixXf bla = MatrixXf::Zero(2, 2);

	bla(0,0) = 13.0;
	bla(0,1) = 5.0;
	bla(1,0) = 2.0;
	bla(1,1) = 4.0;

    EigenSolver<MatrixXf> m_solve(bla);

    VectorXf eigenvalues = VectorXf::Zero(2);
	eigenvalues = m_solve.eigenvalues().real();

	std::cout << "EVs: " << eigenvalues(0) << " and " << eigenvalues(1);

    MatrixXf eigenVectors = MatrixXf::Zero(2, 2);  // matrix (n x m) (points, dims)
    eigenVectors = m_solve.eigenvectors().real();

    tgt::vec2 ev(eigenVectors(0, 0), eigenVectors(0, 1));
    tgt::vec2 ev2(eigenVectors(1, 0), eigenVectors(1, 1));

    std::cout << "\n\n\nEV: ( " << ev.x << " , " << ev.y << " )";
    std::cout << "\nEV2: ( " << ev2.x << " , " << ev2.y << " )\n\n";

    std::cout << "\n\nEigenvalues:\n" << eigenvalues;
    std::cout << "\n\nEigenvectors:\n"  << eigenVectors;
*/



}


Processor* VolumeGenerator::create() const {
    return new VolumeGenerator();
}

/**
 * What to do here? Now we do everything through the update button.
 * We don't want to create new volumes everytime we change any parameter,
 * so just leave this empty and only update when user presses button.
 */
void VolumeGenerator::process() {

    // check for OpenGL errors
    LGL_ERROR;
}


void VolumeGenerator::initialize() throw (VoreenException) {
    // call superclass function first
	VolumeProcessor::initialize();
	volumefactory = new VolumeFactory();
	collection_ = new VolumeCollection();
}

void VolumeGenerator::deinitialize() throw (VoreenException) {

    // call superclass function last
	VolumeProcessor::deinitialize();
}


/**
 * When any of the dimensions is adjusted, check if user wants 'cubed' dimensions.
 * If so, adjust the other two dimensions accordingly.
 */
void VolumeGenerator::onXDimensionsChange(){
	if(dimensions_cubed_prop.get()){
		dimension_y_prop.set(dimension_x_prop.get());
		dimension_z_prop.set(dimension_x_prop.get());
	}
}
void VolumeGenerator::onYDimensionsChange(){
	if(dimensions_cubed_prop.get()){
		dimension_x_prop.set(dimension_y_prop.get());
		dimension_z_prop.set(dimension_y_prop.get());
	}
}
void VolumeGenerator::onZDimensionsChange(){
	if(dimensions_cubed_prop.get()){
		dimension_x_prop.set(dimension_z_prop.get());
		dimension_y_prop.set(dimension_z_prop.get());
	}
}

void VolumeGenerator::update() {
	collection_->clear();

	for(int i=0; i < num_scalar_volumes_prop.get(); i++) {
		if(i % 2 == 0){
			collection_->add(generateScalarVolumeSimple());
		}else{
			collection_->add(generateScalarVolume2());
		}
	}

	for(int i=0; i < num_vector_volumes_prop.get(); i++) {
		if(i % 2 == 0){
			//collection_->add(generateVectorVolumeCircle(0.70, 0.70, 0.05));
			collection_->add(saddleVectorVolume());
		}else{
			collection_->add(generateVectorVolumeCircle(0.75, 0.40, 0.15));
		}
	}

	for(int i=0; i < num_tensor_volumes_prop.get(); i++) {
		collection_->add(generateTensorVolume());
	}
	outport.setData(collection_, true);
}

Volume* VolumeGenerator::generateScalarVolumeSimple(){
	VolumeRAM* scalarVR = volumefactory->create("float", svec3(dimension_x_prop.get(), dimension_y_prop.get(), dimension_z_prop.get()));
	float val = 0.0;
	for(int x=0; x<dimension_x_prop.get(); x++){
		for(int y=0; y<dimension_y_prop.get(); y++){
			for(int z=0; z<dimension_z_prop.get(); z++){
				val = float(y) / float(dimension_y_prop.get());
				scalarVR->setVoxelNormalized(val, x, y, z);
			}
		}
	}
	Volume* result = new Volume(scalarVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}


Volume* VolumeGenerator::generateScalarVolume(){

	VolumeRAM* scalarVR = volumefactory->create("uint8", svec3(dimension_x_prop.get(), dimension_y_prop.get(), dimension_z_prop.get()));
	float val = 0.0;

	// get the shortest dimension and divide by two
	float minimum = std::min(std::min((dimension_x_prop.get()/2), (dimension_y_prop.get()/2)), (dimension_z_prop.get()/2));

	// calculate distances to inner and outter sphere
	float dist_closer = 0.5 * minimum;
	float dist_further = 0.85 * minimum;

	// calculate center of the volume
	float x_center = dimension_x_prop.get() / 2;
	float y_center = dimension_y_prop.get() / 2;
	float z_center = dimension_z_prop.get() / 2;

	// calculate middle between the two spheres
	float cutoff = dist_closer + ((dist_further - dist_closer) / 2);

	float dist_to_center;
	float dist_to_circle;

	// calculate scalar value for each voxel in volume
	for(int x=0; x<dimension_x_prop.get(); x++){
		for(int y=0; y<dimension_y_prop.get(); y++){
			for(int z=0; z<dimension_z_prop.get(); z++){

				dist_to_center = sqrt( (x - x_center)*(x - x_center) + (y - y_center)*(y - y_center) + (z - z_center)*(z - z_center) );

				if(dist_to_center <= cutoff){
					// first (smaller) peak: max .5*255: 0.5 / (1 + x^2)
					dist_to_circle = abs(dist_closer - dist_to_center);
					val = 0.5 / (1.0 + (dist_to_circle * dist_to_circle));
				}else{
					// second (larger) peak: max 255: 1.0 / (1 + x^2)
					dist_to_circle = abs(dist_further - dist_to_center);
					val = 1.0 / (1.0 + (dist_to_circle * dist_to_circle));

					// small hack to ensure inner circle has a 'unique' value between 0.35 and 0.50
					if(val <= 0.51 && val >= 0.35){
						if(val <= 0.43){ val = 0.35; }else{ val = 0.51; }
					}
				}

				scalarVR->setVoxelNormalized(val, x, y, z);
			}
		}
	}

	Volume* result = new Volume(scalarVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}

Volume* VolumeGenerator::generateScalarVolume2(){

	VolumeRAM* scalarVR = volumefactory->create("uint8", svec3(dimension_x_prop.get(), dimension_y_prop.get(), dimension_z_prop.get()));
	float val = 0.0;

	// get the shortest dimension and divide by two
	float minimum = std::min(std::min((dimension_x_prop.get()/2), (dimension_y_prop.get()/2)), (dimension_z_prop.get()/2));

	// calculate center of the volume
	float x_14 = dimension_x_prop.get() * 0.14;
	float x_16 = dimension_x_prop.get() * 0.16;
	float x_84 = dimension_x_prop.get() * 0.84;
	float x_86 = dimension_x_prop.get() * 0.86;
	float y_47 = dimension_y_prop.get() * 0.47;
	float y_53 = dimension_y_prop.get() * 0.53;

	// calculate scalar value for each voxel in volume
	for(int x=0; x<dimension_x_prop.get(); x++){
		for(int y=0; y<dimension_y_prop.get(); y++){
			for(int z=0; z<dimension_z_prop.get(); z++){

				if(y > y_47 && y < y_53){
					 val = 0.9;
				}else if( (x > x_14 && x < x_16) || (x > x_84 && x < x_86) ){
					val = 0.5;
				}else{
					val = 0.0;
				}

				scalarVR->setVoxelNormalized(val, x, y, z);
			}
		}
	}

	Volume* result = new Volume(scalarVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}

Volume* VolumeGenerator::generateScalarVolumeHelix(){

	VolumeRAM* scalarVR = volumefactory->create("uint8", svec3(dimension_x_prop.get(), dimension_y_prop.get(), dimension_z_prop.get()));
	float val = 0.0;

	// calculate scalar value for each voxel in volume
	for(int x=0; x<dimension_x_prop.get(); x++){
		for(int y=0; y<dimension_y_prop.get(); y++){
			for(int z=0; z<dimension_z_prop.get(); z++){
				val = 1 / ( (x - cos((float)z))*(x - cos((float)z)) + (y - sin((float)z))*(y - sin((float)z)) );

				scalarVR->setVoxelNormalized(val, x, y, z);
			}
		}
	}

	Volume* result = new Volume(scalarVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}

Volume* VolumeGenerator::generateVectorVolume(){
	/**
	 * TODO
	 * 	- Write a more interesting vector field (look up Hairy Ball Theorem on wikipedia)
	 * 	- Look at FlowRenderer2D / FlowRenderer3D for example implementations of Spot Noise, LIC etc
	 * 	- Look at GeometryProcessor to see how to do an input for multiple sources (so I can get an array of Volumes and FlowreenAdapter output things as input to my processor)
	 */

	// NB: for this data to be compatible with flowreen: use "Vector3(float)" and run it through a FlowreenAdapter module!
	VolumeRAM* vectorVR = volumefactory->create("Vector3(float)", svec3(dimension_x_prop.get(), dimension_y_prop.get(), dimension_z_prop.get()));
	float valx = 0.0;
	float valy = 0.0;
	float valz = 0.0;

	float minimum = std::min(std::min((dimension_x_prop.get()/2), (dimension_y_prop.get()/2)), (dimension_z_prop.get()/2));

	float dist_closer = 0.5 * minimum;
	float dist_further = 0.85 * minimum;

	float vec_to_center_x = 0.0;
	float vec_to_center_y = 0.0;
	float vec_to_center_z = 0.0;

	float x_center = dimension_x_prop.get() / 2;
	float y_center = dimension_y_prop.get() / 2;
	float z_center = dimension_z_prop.get() / 2;

	float cutoff = dist_closer + ((dist_further - dist_closer) / 2);

	float dist_to_center;
	float dist_to_circle;

	for(int x=0; x<dimension_x_prop.get(); x++){
		for(int y=0; y<dimension_y_prop.get(); y++){
			for(int z=0; z<dimension_z_prop.get(); z++){

				vec_to_center_x = x_center - x;
				vec_to_center_y = y_center - y;
				vec_to_center_z = z_center - z;

				dist_to_center = sqrt( (vec_to_center_x)*(vec_to_center_x) + (vec_to_center_y)*(vec_to_center_y) + (vec_to_center_z)*(vec_to_center_z) );

				vec_to_center_x = vec_to_center_x / dist_to_center;
				vec_to_center_y = vec_to_center_y / dist_to_center;
				vec_to_center_z = vec_to_center_z / dist_to_center;

				/*
				 * See tgt's matrix.h for multiplication of vector aiming to center with a rotation matrix4
				 * (should it be a matrix3?):
				 *
/// This operator is not available in GLSL but very useful: A mat4 * vec3 operator, returning a vec3
/// note: the resulting vec3 is divided by the w component of the temporary result!
template<class T>
Vector3<T> operator * (const Matrix4<T>& m, const Vector3<T>& v) {
    Vector4<T> v4(v, 1);
    v4 = m * v4;
    v4 /= v4.w;
    return Vector3<T>(v4.elem);
}
				 *
				 */

/*
				if(dist_to_center <= cutoff){
					// first (smaller) peak: max .5*255: 0.5 / (1 + x^2)
					dist_to_circle = abs(dist_closer - dist_to_center);
					val = 0.5 / (1.0 + (dist_to_circle * dist_to_circle));
				}else{
					// second (larger) peak: max 255: 1.0 / (1 + x^2)
					dist_to_circle = abs(dist_further - dist_to_center);
					val = 1.0 / (1.0 + (dist_to_circle * dist_to_circle));

					// small hack to ensure inner circle has a 'unique' value between 0.35 and 0.50
					if(val <= 0.51 && val >= 0.35){
						if(val <= 0.43){ val = 0.35; }else{ val = 0.51; }
					}
				}
*/
				vectorVR->setVoxelNormalized(vec_to_center_x, x, y, z, 0);
				vectorVR->setVoxelNormalized(vec_to_center_y, x, y, z, 1);
				vectorVR->setVoxelNormalized(vec_to_center_z, x, y, z, 2);
			}
		}
	}

	Volume* result = new Volume(vectorVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
    // Volume(VolumeRepresentation* const volume, const tgt::vec3& spacing, const tgt::vec3& offset, const tgt::mat4& transformation = tgt::mat4::identity);
	return result;
}

Volume* VolumeGenerator::generateVectorVolume2(){
	// TODO generate a vector volume that corresponds with generateScalarVolume2()
	VolumeRAM* vectorVR = volumefactory->create("Vector3(float)", svec3(dimension_x_prop.get(), dimension_y_prop.get(), dimension_z_prop.get()));
	float valx = 0.0;
	float valy = 0.0;
	float valz = 0.0;

	float x_14 = dimension_x_prop.get() * 0.14;
	float x_16 = dimension_x_prop.get() * 0.16;
	float x_84 = dimension_x_prop.get() * 0.84;
	float x_86 = dimension_x_prop.get() * 0.86;
	float y_47 = dimension_y_prop.get() * 0.47;
	float y_53 = dimension_y_prop.get() * 0.53;

	for(int x=0; x<dimension_x_prop.get(); x++){
		for(int y=0; y<dimension_y_prop.get(); y++){
			for(int z=0; z<dimension_z_prop.get(); z++){


				if(y == floor(y_47)){
					valx = -1.0;
					valy = 0.0;
					valz = 0.0;
				}else if( y == ceil(y_53) ){
					valx = 1.0;
					valy = 0.0;
					valz = 0.0;
				}else if( x == floor(x_14) || x == floor(x_84)){
					valx = 0.0;
					valy = 1.0;
					valz = 0.0;
				}else if( x == ceil(x_16) || x == ceil(x_86)){
					valx = 0.0;
					valy = -1.0;
					valz = 0.0;
				}else{
					valx = 0.0;
					valy = 0.0;
					valz = 0.0;
				}

				vectorVR->setVoxelNormalized(valx, x, y, z, 0);
				vectorVR->setVoxelNormalized(valy, x, y, z, 1);
				vectorVR->setVoxelNormalized(valz, x, y, z, 2);
			}
		}
	}

	Volume* result = new Volume(vectorVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}

/**
 * Returns a vector volume with only straight vectors in one direction
 */
Volume* VolumeGenerator::generateVectorVolumeWaves(){
	VolumeRAM* vectorVR = volumefactory->create("Vector3(float)", svec3(dimension_x_prop.get(), dimension_y_prop.get(), dimension_z_prop.get()));
	float valx = 0.0;
	float valy = 1.0;
	float valz = 0.0;
	float angular_velocity = 0.15;

	tgt::vec3 singularity(ceil(0.40*dimension_x_prop.get()), ceil(0.70*dimension_y_prop.get()), ceil(0.50*dimension_z_prop.get()));

	for(int x=0; x<dimension_x_prop.get(); x++){
		for(int y=0; y<dimension_y_prop.get(); y++){
			for(int z=0; z<dimension_z_prop.get(); z++){

				//FIXME this seems to just get a value of -0.2, LIC shows lines directed to top-left

				float dist = calculateVectorMagnitude(tgt::vec3(singularity.x - x, singularity.y - y, 0.0));

				if(dist > 20){
					valx = sin((float)y / (float)5.0);
				}else{
					valx = (-1.0 * angular_velocity * (y - singularity.y));
					valy = dist / (float)20;
				}

				vectorVR->setVoxelNormalized(valx, x, y, z, 0);
				vectorVR->setVoxelNormalized(valy, x, y, z, 1);
				vectorVR->setVoxelNormalized(valz, x, y, z, 2);
			}
		}
	}

	Volume* result = new Volume(vectorVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}

/**
 * Returns a vector volume with some circles. Vector volume is identical in 1 direction,
 * i.e. in the z-dimension, every 2D slice looks the same. This means the z-component is always 0.
 */
Volume* VolumeGenerator::generateVectorVolumeCircles(){
	VolumeRAM* vectorVR = volumefactory->create("Vector3(float)", svec3(dimension_x_prop.get(), dimension_y_prop.get(), dimension_z_prop.get()));
	float valx = 0.0;
	float valy = 0.0;
	float valz = 0.0;

	float angular_velocity = 0.15;
	float pointA_x = 0.25*dimension_x_prop.get();
	float pointA_y = 0.75*dimension_y_prop.get();
	float pointB_x = 0.75*dimension_x_prop.get();
	float pointB_y = 0.25*dimension_y_prop.get();

	for(int x=0; x<dimension_x_prop.get(); x++){
		for(int y=0; y<dimension_y_prop.get(); y++){
			for(int z=0; z<dimension_z_prop.get(); z++){
				//TODO calculate values here

				if(y > x){ // rotational flow around point A
					valx = -1.0 * angular_velocity * (y - pointA_y);
					valy = angular_velocity * (x - pointA_x);
				}else if(x > y){ // rotational flow around point B
					valx = -1.0 * angular_velocity * (y - pointB_y);
					valy = angular_velocity * (x - pointB_x);
				}else{ // x == y, just a straight vector along diagonal.
					valx = 1.0;
					valy = 1.0;
				}

				vectorVR->setVoxelNormalized(valx, x, y, z, 0);
				vectorVR->setVoxelNormalized(valy, x, y, z, 1);
				vectorVR->setVoxelNormalized(valz, x, y, z, 2);
			}
		}
	}

	Volume* result = new Volume(vectorVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}

Volume* VolumeGenerator::generateVectorVolumeCircle(float ax, float ay, float angular_velocity){
	VolumeRAM* vectorVR = volumefactory->create("Vector3(float)", svec3(dimension_x_prop.get(), dimension_y_prop.get(), dimension_z_prop.get()));
	float valx = 0.0;
	float valy = 0.0;
	float valz = 0.0;

	float pointA_x = ax*dimension_x_prop.get();
	float pointA_y = ay*dimension_y_prop.get();

	for(int x=0; x<dimension_x_prop.get(); x++){
		for(int y=0; y<dimension_y_prop.get(); y++){
			for(int z=0; z<dimension_z_prop.get(); z++){
				//TODO perhaps invert magnitudes (ie larger values close to nodes)
				valx = (-1.0 * angular_velocity * (y - pointA_y));
				valy = (angular_velocity * (x - pointA_x));

				vectorVR->setVoxelNormalized(valx, x, y, z, 0);
				vectorVR->setVoxelNormalized(valy, x, y, z, 1);
				vectorVR->setVoxelNormalized(valz, x, y, z, 2);
			}
		}
	}

	Volume* result = new Volume(vectorVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}

Volume* VolumeGenerator::generateVectorVolumeSingularities(){
	VolumeRAM* vectorVR = volumefactory->create("Vector3(float)", svec3(dimension_x_prop.get(), dimension_y_prop.get(), dimension_z_prop.get()));
	float ax = 0.0;
	float ay = 0.0;
	float az = 0.0;

	float bx = 0.0;
	float by = 0.0;
	float bz = 0.0;

	float crossx = 0.0;
	float crossy = 0.0;
	float crossz = 0.0;

	float angular_velocity = 0.15;
	float pointA_x = ceil(0.40*dimension_x_prop.get());
	float pointA_y = ceil(0.70*dimension_y_prop.get());
	float pointA_z = 5.0;
	float pointB_x = ceil(0.75*dimension_x_prop.get());
	float pointB_y = ceil(0.25*dimension_y_prop.get());
	float pointB_z = 10.0;

	for(int x=0; x<dimension_x_prop.get(); x++){
		for(int y=0; y<dimension_y_prop.get(); y++){
			for(int z=0; z<dimension_z_prop.get(); z++){


				/**
				 * TODO
				 * - Figure out how to define a 3D vector field with singularities!
				 * - Once I've got that figured out, define two of them and use the cross-product (I think)
				 * 	 of those, to get a vector-field with 2 singularities.
				 * - Look into Saddle points, Repelling/Attracting nodes, Repelling and Attracting focus,
				 * 	 and Center (Topological Representations of Vector Fields, paper by Theisel et al.)
				 */

				// currently: cross-product of repelling node and center and attracting nodes
				ax = (x - pointA_x) * 0.1;
				ay = (y - pointA_y) * 0.1;
				az = (z - pointA_z) * 0.1;

				bx = (x - pointB_x) * 0.1;
				by = (y - pointB_y) * 0.1;
				bz = (z - pointB_z) * 0.1;

				//<a,b,c> x <d,e,f> = <bf-ce, cd-af, ae-bd>
				crossx = ay*bz - az*by;
				crossy = az*bx - ax*bz;
				crossz = ax*by - ay*bx;

				vectorVR->setVoxelNormalized(crossx, x, y, z, 0);
				vectorVR->setVoxelNormalized(crossy, x, y, z, 1);
				vectorVR->setVoxelNormalized(crossz, x, y, z, 2);
			}
		}
	}

	Volume* result = new Volume(vectorVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}


Volume* VolumeGenerator::generateVectorVolume100(){
	VolumeRAM* vectorVR = volumefactory->create("Vector3(float)", svec3(dimension_x_prop.get(), dimension_y_prop.get(), dimension_z_prop.get()));

	tgt::vec3 vect(0.0, 0.0, 0.0);

	tgt::vec3 singularity(ceil(0.40*dimension_x_prop.get()), ceil(0.70*dimension_y_prop.get()), ceil(0.50*dimension_z_prop.get()));

	float max_distance = calculateVectorMagnitude( tgt::vec3(
			(dimension_x_prop.get() - singularity.x),
			singularity.y,
			singularity.z ) );
	float dist_factor = 0.0;

	for(int x=0; x<dimension_x_prop.get(); x++){
		for(int y=0; y<dimension_y_prop.get(); y++){
			for(int z=0; z<dimension_z_prop.get(); z++){

				vect.x = singularity.x - x;
				vect.y = singularity.y - y;
				vect.z = singularity.z - z;

				dist_factor = calculateVectorMagnitude(vect);
				if(dist_factor < 0.0001){
					dist_factor = 0.0001;
				}
				dist_factor = dist_factor / (float)10;


				vect = normalizeVector(vect);

				vect.x = vect.x / (dist_factor);
				vect.y = vect.y / (dist_factor);
				vect.z = vect.z / (dist_factor);

				vectorVR->setVoxelNormalized(vect.x, x, y, z, 0);
				vectorVR->setVoxelNormalized(vect.y, x, y, z, 1);
				vectorVR->setVoxelNormalized(vect.z, x, y, z, 2);
			}
		}
	}

	Volume* result = new Volume(vectorVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}

Volume* VolumeGenerator::generateVectorVolume200(){
	VolumeRAM* vectorVR = volumefactory->create("Vector3(float)", svec3(dimension_x_prop.get(), dimension_y_prop.get(), dimension_z_prop.get()));

	tgt::vec3 vect(0.0, 0.0, 0.0);

	tgt::vec3 singularity(ceil(0.80*dimension_x_prop.get()), ceil(0.20*dimension_y_prop.get()), ceil(0.50*dimension_z_prop.get()));

	float max_distance = calculateVectorMagnitude(
			tgt::vec3( singularity.x,
					   (dimension_y_prop.get() - singularity.y),
					   singularity.z ) );
	float dist_factor = 0.0;

	for(int x=0; x<dimension_x_prop.get(); x++){
		for(int y=0; y<dimension_y_prop.get(); y++){
			for(int z=0; z<dimension_z_prop.get(); z++){

				vect.x = singularity.x - x;
				vect.y = singularity.y - y;
				vect.z = singularity.z - z;

				dist_factor = calculateVectorMagnitude(vect);
				if(dist_factor < 0.1){
					dist_factor = 0.1;
				}
				dist_factor = dist_factor / (float)10;

				vect = normalizeVector(vect);

				if(dist_factor != 0.0){
					vect.x = vect.x / (dist_factor);
					vect.y = vect.y / (dist_factor);
					vect.z = vect.z / (dist_factor);
				}

				vectorVR->setVoxelNormalized(vect.x, x, y, z, 0);
				vectorVR->setVoxelNormalized(vect.y, x, y, z, 1);
				vectorVR->setVoxelNormalized(vect.z, x, y, z, 2);

			}
		}
	}

	Volume* result = new Volume(vectorVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}

Volume* VolumeGenerator::saddleVectorVolume(){

	int xdim = dimension_x_prop.get();
	int ydim = dimension_y_prop.get();
	int zdim = dimension_z_prop.get();

	// 1. Calculate Tensor field with saddle
	float v11 = 1.0;
	float v12 = 0.0;
	float v21 = 0.0;
	float v22 = 1.0;

	tgt::ivec2 saddle(int(xdim*0.10),int(ydim*0.30));
	float dist = 0.0;

	//TODO is this meant to be between 0 and 1?
	float decay = 0.5;
	float factor = 0.0;
/*
	CTensorField *tfield = new CTensorField(xdim, ydim);
	tfield->init();
	for(int x=0; x<xdim; x++){
		for(int y=0; y<ydim; y++){

			tgt::ivec2 v(saddle.x-x, saddle.y-y);
			dist = sqrtf(tgt::dot(v, v));

			factor = exp(-1 * decay * (dist*dist));

			//TODO calculate v11 v12 v21 v22
			v11 = factor * (x*x - y*y);
			v12 = factor * (-2*x*y);
			v21 = v12;
			v22 = -v11;

			tfield->set(x, y, new CTensor(v11, v12, v21, v22));
		}
	}
*/

	// 2. calculate eigenvectors per tensor
	VolumeRAM* vectorVR = volumefactory->create("Vector3(float)", svec3(xdim, ydim, zdim));
	tgt::vec2 vect(0.0, 0.0);
	for(int x=0; x<xdim; x++){
		for(int y=0; y<ydim; y++){
			//vect = tfield->get(x,y).eigenVector3D();
			tgt::ivec2 v(saddle.x-x, saddle.y-y);
			dist = sqrtf(tgt::dot(v, v));
			factor = exp(-1 * decay * (dist*dist));

			v11 = factor * (x*x - y*y);
			v12 = factor * (-2*x*y);
			v21 = v12;
			v22 = -v11;

			bool show = x==y && x < int(xdim*0.10);

			vect = calcEigenvector(tgt::vec2(v11, v21), tgt::vec2(v12, v22), show);
			for(int z=0; z<zdim; z++){
				vectorVR->setVoxelNormalized(vect.x, x, y, z, 0);
				vectorVR->setVoxelNormalized(vect.y, x, y, z, 1);
				vectorVR->setVoxelNormalized(0.0, x, y, z, 2);
			}
		}
	}

	Volume* result = new Volume(vectorVR, vec3(1.0, 1.0, 1.0), vec3(0.0, 0.0, 0.0));
	centerVolume(result);
	return result;
}

Volume* VolumeGenerator::generateTensorVolume(){
	// TODO
	return generateScalarVolume();
}

float VolumeGenerator::calculateVectorMagnitude( tgt::vec3 v ) {
	return sqrtf(tgt::dot(v, v));
}

tgt::vec3 VolumeGenerator::normalizeVector( tgt::vec3 v ) {
	return v / calculateVectorMagnitude(v);
}

/**
 * Helper function to calculate eigenvector of 2x2 matrix with
 * c1 and c2 as first and second column respectively
 */
tgt::vec2 VolumeGenerator::calcEigenvector(tgt::vec2 c1, tgt::vec2 c2, bool show){

	using namespace Eigen;

	MatrixXf tmp = MatrixXf::Zero(2, 2);

	tmp(0,0) = c1.x;
	tmp(1,0) = c1.y;
	tmp(0,1) = c2.x;
	tmp(1,1) = c2.y;

    EigenSolver<MatrixXf> solver(tmp);

    MatrixXf eigenVectors = MatrixXf::Zero(2, 2);
    eigenVectors = solver.eigenvectors().real();


    tgt::vec2 ev(eigenVectors(1, 0), eigenVectors(1, 1));

    if(!(ev.x == 0 && ev.y == 1) || show){
    	std::cout << "\nMatrix: " << tmp << " - vectors: " << eigenVectors << " - ev: " << ev;
    }
    return ev;
}

} // namespace
